Calculations of the Structure of Basin Volumes for Mechanically Stable Packings 



O 

(N 

O 

Q 

On 



O 

a 

C 

o 
o 



> 

(N 



X 



S. S. Ashwin 1 ' 2 , Jerzy Blawzdziewicz 3 , Corey S. O'Hern 1 ' 2 , and Mark D. Shattuck 4 
1 Department of Mechanical Engineering and Materials Science, Yale University, New Haven, CT 06520-8286 
2 Department of Physics, Yale University, New Haven, CT 06520-8120 
3 Department of Mechanical Engineering, Texas Tech University, Lubbock, TX 74-909-1021 and 
4 Benjamin Levich Institute and Physics Department, 
The City College of the City University of New York, New York, NY 10031 
(Dated: December 20, 2011) 

There are a finite number of distinct mechanically stable (MS) packings in model granular systems 
composed of frictionless spherical grains. For typical packing-generation protocols employed in 
experimental and numerical studies, the probabilities with which the MS packings occur are highly 
nonuniform and depend strongly on parameters in the protocol. Despite intense work, it is extremely 
difficult to predict a priori the MS packing probabilities, or even which MS packings will be the most 
versus the least probable. We describe a novel computational method for calculating the MS packing 
probabilities by directly measuring the volume of the MS packing 'basin of attraction', which we 
define as the collection of initial points in configuration space at zero packing fraction that map to a 
given MS packing by following a particular dynamics in the density landscape. We show that there is 
a small core region with volume V„ surrounding each MS packing n in configuration space in which all 
initial conditions map to a given MS packing. However, we find that the MS packing probabilities are 
very weakly correlated with core volumes. Instead, MS packing probabilities obtained using initially 
dilute configurations are determined by complex geometric features of the basin of attraction that 
are distant from the MS packing. 

PACS numbers: 63.50.Lm,83.80.Fg61.43.-j,64.70.ps 



I. INTRODUCTION 

In contrast to equilibrium, thermal systems, the struc- 
tural and mechanical properties of dense granular ma- 
terials and other athermal particulate systems depend 
strongly on the protocol used to create them. For ex- 
ample, a number of studies have shown that the packing 
fraction of granular assemblies can vary from values asso- 
ciated with random loose [l[ to random close packing Q 
as a function of the vibration amplitude and tapping his- 
tory |H, HJ. In addition, the force chain networks that 
form, and thus the shear modulus of granular packings 
depend on whether they have been generated via shear, 
isotropic compression [5|, or sedimentation via single- 
particle or collective deposition @. 

The protocol dependence in dense granular systems 
arises from the nonlinear, dissipative, and frictional con- 
tact interactions between grains @- Despite active re- 
search in this area, the distinct contributions from each 
of these interactions to protocol dependence has not been 
determined. In this manuscript, we will investigate the 
protocol-dependence of static granular packings by fo- 
cusing on a simple system of frictionless spherical parti- 
cles that interact via purely repulsive linear spring and 
velocity-dependent damping forces. For a fixed set of 
boundary conditions, there are finite number of distinct 
mechanically stable (MS) packings of frictionless parti- 
cles, which grows exponentially with the number of parti- 
cles N [U . MS packings exist as discrete points in config- 
uration space that are characterized by the packing frac- 
tion <j)j and N particle coordinates Rj = {fx, r*2 . . . r/v} 
and coincide with local minima of the density land- 
scape 0, H3 (or local minima of the potential energy 



landscape with zero potential V = 0). We have shown 
recently in both simulations and experiments that the 
probabilities with which these distinct MS packings oc- 
cur are highly nonuniform and depend on parameters of 
the packing-generation protocol including the compres- 
sion rate, damping coefficient, and initial packing frac- 
tion |ITl[l2j. However, one cannot yet determine a priori 
which MS packings are the most versus the least prob- 
able, much less calculate the packing probabilities as a 
function of the packing-generation protocol. 

Here we describe a novel computational method for 
calculating the MS packing probabilities by directly mea- 
suring the volume of the MS packing 'basin of attrac- 
tion', which we define as the collection of initial points 
in configuration space (i.e., the red region in Fig. [T]) at 
zero packing fraction that map to a given MS packing 
by following a particular dynamics in the density land- 
scape. Note that our definition of the basin of attraction 
is protocol-dependent, and thus the basin volume will 
vary with the rate at which energy is dissipated, the com- 
pression rate, and other parameters. In contrast, basins 
of attraction for glassy systems [l3|, [13] are defined as the 
set of initial dense liquid configurations that map to the 
'nearest' local minimum using steepest descent dynamics 
at fixed density. Our definition of basin volumes is more 
relevant for granular systems, in which MS packings are 
generated from initially dilute configurations. 

To aid in the calculation of the basin volumes, we intro- 
duce the unweighted basin profile function /£(/), which is 
the fraction of points on a hypersurface in configuration 
space a distance I from the nth MS packing that maps via 
a given dynamics (labeled T) to MS packing n. We will 
show that there is a hyperspherical core region surround- 
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FIG. 1: (color online) (top) A schematic of the basin of attrac- 
tion (red) in diV-dimensional configuration space for a typi- 
cal mechanically stable (MS) packing (black dot), (bottom) 
The corresponding unweighted basin profile function f^(l) is 
plotted as a function of distance I from MS packing n for 
packing-generation protocol T. f„(l) begins to decay from 1 
beyond an approximately spherically symmetric core size l c , 
while for / > l c the basin is highly branched, thread-like, and 

/£(0->o- 

ing each MS packing in which f^(l) = 1 for I < l c , while 
further from the MS packing, the basin becomes highly 
branched, thread-like, and f F (l) 0. (See Fig.[TJ) This 
picture raises several key questions: 1) Are the MS pack- 
ing packing probabilities determined by the size l c of the 
core region in configuration space or dominated by contri- 
butions from the thread-like regions, and 2) do the mor- 
phologies of the basins of attraction depend sensitively 
on protocol? We will show below that the MS packing 
probabilities are not strongly correlated with the volume 
of the core regions in configuration space and are instead 
determined by features of the density landscape that are 
far from each MS packing packing. Thus, novel compu- 
tational geometry techniques dH must be developed to 
understand the key features of configuration space that 
control MS packing probabilities. 



II. METHODS 

To perform our calculations of basin volumes, we fo- 
cused on a well-characterized model system composed of 
N frictionless disks in 2D that interact via purely repul- 
sive linear spring and velocity-dependent damping forces. 
N is varied from 3 to 100, and the particles are enclosed 
in a square cell with fixed walls of length L = 1. In- 
teractions with the walls match those between the parti- 
cles. We consider both monodisperse and bidisperse sys- 
tems, where the bidisperse mixtures contain half large 
and half small disks (N s — Ni — N/2) with diameter 
ratio ai/a s = 1.4. 



(b) 



(c) 



FIG. 2: (color online) The microstate basins of attraction for 
a system of three monodisperse frictionless disks, where par- 
ticles 2 and 3 are initially located at positions (0.2, 0.6) and 
(0.45, 0.85) in the x-y plane (with the origin in the lower left 
corner). Results are shown for two damping coefficients, (a) 
6=1 and (b) b = 0.1. The position of each pixel represents 
the initial position of particle 1 and its color corresponds to 
one of the 11 out of 24 microstates in (c) to which the system 
evolved under the compression protocol. For N = 3 monodis- 
perse systems, there is one distinct MS packing (Ms = 1) with 
four polarizations (hue; rows) and six permutations (satura- 
tion; columns) for a total of Mm = 24 microstates. 

In a number of previous studies, we described the 
MS 'packing finder' that generates a mechanically stable 
packing via isotropic compression at <f>j with infinitesimal 
overlap from an arbitrary initial condition at (f> = [ill ]. 
Briefly, the algorithm includes the following steps. For 
each trial, we initialize the system with random particle 
positions inside the unit square at 4> — and zero veloci- 
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TABLE I: The number of distinct mechanically stable pack- 
ings Ma and total number of microstates M m versus the num- 
ber of particles N. For JV = 3 we consider monodisperse 
systems. For the other system sizes, results are given for 
bidisperse mixtures. For N — 12 we estimate M a and Mm- 
We do not include unstable packings such as the one in the 
upper left corner of Fig. [3] in which the 'rigid backbone' of 
particles can translate. 

ties. We then compress the system in steps of Acj> = 10~ 4 
and relax the small particle overlaps after each step by 
solving Newton's equations of motion with damping, 



= ^F(r i3 ) - bv, 



(1) 



where m, er, and a, are the particle mass, diameter, and 
acceleration, 



(2) 



e is the characteristic energy of the repulsive spring inter- 
action, 0(x) is the Heaviside step function, b = ba/^/me 
is the damping coefficient, fjj is the unit vector con- 
necting the centers of particles i and j and r»j is their 
separation, until the kinetic energy per particle falls be- 
low a specified tolerance K/eN < K to \ = 10~ 25 . We 
studied a wide range of values for the damping coeffi- 
cient from b = 10 -2 to 10, which mimics steepest descent 
dynamics. The packing-generation algorithm terminates 
when the minimized total potential energy per particle 
V/eN > Vtoi = 10~ 16 . As in previous studies on similar 
systems with periodic boundary conditions, we distin- 
guish MS packings based on the spectrum of nontrivial 
eigenvalues of the dynamical matrix [l2j . and we hnd 
that the number of distinct MS packings M s grows expo- 
nentially with TV as shown in Table fl] The M s = 6 and 
80 distinct MS packings for N — 4 and 6 are shown in 
Figs. [3] and |U The packing finder does produce a small 
number of unstable packings as shown in Fig. 02 but these 
are not included in the analyses. 

The fundamental quantity in our approach is the un- 
weighted basin profile function defined as 



f n {l) = / dRG r (R, R n j)S(\R - R n 7 \ - I), 



(3) 



where is sampled on hyperspherical shells a distance 
I from MS packing n, T is the specified compression dy- 
namics, 8(x) is the Dirac delta function, Gr(R, Rj) = 1 
for points R in configuration space that map to MS pack- 
ing R 7 j, and otherwise. As an illustrative example, 




FIG. 3: (color online) The bottom six configurations are the 
Ms = 6 distinct mechanically stable packings for bidisperse 
systems with N = 4. (We will refer to configurations 1 
through 6 counting in ascending order from left to right and 
bottom to top.) The particles shaded blue form the force- 
bearing backbone of the mechanically stable packing. Parti- 
cles shaded red are 'rattlers' with fewer than 3 contacts. The 
packing finder generates a small number of unstable config- 
urations similar to that shown in the upper left corner with 
probability less than 0.2%, but these are not included in the 
analyses. 



we calculate slices of Gr{R,Rj) for N = 3, which has 
a single MS packing with M m microstates — 6 particle- 
label permutations and 4 polarizations obtained by ap- 
plying all possible reflections and rotations in 2D consis- 
tent with the square cell boundary conditions [l6j |. In 
Fig. [21 we plot the microstate basins of the attraction 

Y^=i n G{{ri^2,^},R n j) for fixed r% = (0.2,0.6) and 
f% = (0.45,0.85). 

We calculate the unweighted basin profile function 
fn{l) using two procedures; the first method is efficient 
and accurate for small I and the second for large I. For 
method 1, we generate at least M = 10 6 points ran- 
domly on the surface of a 2A r -dimensional hypersphere 
centered on the MS packing with radius /. We then in- 
put each of these configurations as initial configurations 




FIG. 4: (color online) The J\f s — 80 distinct mechanically sta- 
ble packings for bidisperse systems with N = 6. The particles 
shaded blue form the force-bearing backbone of the mechan- 
ically stable packing. Particles shaded red are 'rattlers' with 
fewer than 3 contacts. The unweighted and weighted basin 
profile functions are shown in Fig. [11] for the configuration in 
the third row that is shaded gray. 

into the MS packing finder with packing fraction fa = 0. 
If a given initial condition belongs to the basin of attrac- 
tion of MS packing n, the packing finder will generate 
packing n. Otherwise, the initial condition belongs to 
a different basin. For the system sizes where we can 
achieve complete enumeration, we found that the crite- 
rion, maxi(c^ — di)/di < 10~ 6 , was sufficiently sensitive 
to distinguish MS packings, where d\ is the ith sorted 
eigenvalue of the dynamical matrix for MS packing j. 
From method 1, the unweighted basin profile function 
for MS packing n is 



where M n is the number of initial conditions at I that 
map to to packing n. 

We define the basin volume for MS packing n generated 
using compression dynamics T as 

V n = / S^(l)dl, (5) 
Jo 

where 

Sl(l) = A 2N f^(l)l 2N - 1 V n N s \N l \ (6) 

is the (angle-averaged) weighted basin profile function, 
A k = 27r fc / 2 /r(fc/2) is the surface area of a fc-dimensional 
unit sphere, and V n is the number of distinct polariza- 
tions for MS packing n [l6|. The probability of MS 
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FIG. 5: (color online) [top] The unweighted /n(0 and [bot- 
tom] weighted basin profile functions measured using 
methods 1 (circles) and 2 (solid lines) for MS packings 1 (high- 
est probability; dark blue line) and 4 (lowest probability; light 
red line) shown in Fig. [3] plotted on a log- log scale for N = 4 
and 6=1. The vertical lines indicate l c for each MS packing. 

packing n for a given compression protocol T is pro- 
portional to its basin volume, = /Vtot, where 

Vtot = T * lV T = L ™ = 1 . 

Method 1 becomes extremely inefficient at calculating 
fn(l) f° r large I > l c . Thus, in this regime we implement 
method 2, which was previously employed to calculate 
the probabilities directly Q. For this method, we gen- 
erate at least 10 6 random points in configuration space 
and input these into the packing finder with fa = 0. The 
fraction of random initial configurations that map to MS 
packing n determines P£. We can then calculate 
from P^ using Eqs.[5]and[ni Note that an advantage of 
method 2 is that each initial condition provides informa- 
tion about P^ for some n and for N S \N[\ distances I by 
permuting the labels of the final MS packing. 



III. RESULTS 

Typical basin profile functions /„(/) are shown for the 
most and least probable MS packings (1 and 4 in Fig. [3]) 
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FIG. 6: The volume V„ of the hyperspherical cores surround- 
ing each MS packing n (relative to Vtot) plotted as a function 
of the MS packing probability P„ for each MS packing for 
N = 4 (asterisks) and 6 (circles) obtained using method 1 
with damping parameter b = 1. The solid (dashed) line has 
slope 1 (7.5). 



for N = 4 in the top panel of Fig. [5] For small dis- 
tances from the MS packing I < l c , f n (l) = 1. Beyond 
the core size l c , which can vary strongly from one MS 
packing to another, f„(l) decays rapidly to zero. In the 
bottom panel of Fig. [5J we show the weighted basin pro- 
file S%(1) for the same N = 4 MS packings. Since S£(Z) 
is obtained by multiplying /„(£) by l 2N ~ x , the probabil- 
ities for obtaining MS packings (when starting from zero 
packing fraction) are determined by distances I > l c - For 
N = 4 the average core size is (l c ) » 0.1, the small parti- 
cle diameter is a = 0.3, but the average length scale that 
yields 50% of the packing probabilities (near the peak in 
Sn(l)) is (l p ) 0.5. We will show below that l p grows 
with increasing system size. We have validated the re- 
sults by ensuring that methods 1 and 2 yield the same 
values for f„(l) and S£(Z) over the range in I in which 
the calculations overlap. 

In the top panel of Fig. [SJ we show that the core size for 
the most probable N = 4 MS packing is larger than that 
for the least probable MS packing, which may suggest 
that there is a correlation between the core size and the 
MS packing probabilities. To investigate to what extent 
the hyperspherical core surrounding each MS packing de- 
termines the packing probabilities, we approximate the 
basin volume by the volume of a hypersphere of radius l c , 
V£ = tt n II n /T(N + 1), for each MS packing. In Fig. |1 
we plot V^/Vtot versus P^ for N — 4 and 6. We find 
two key results: 1) The volumes V£/Vt t ar e smaller by 
many orders of magnitude than the probabilities P^ and 
2) A fit to the data for TV = 6 yields V£/V tot ~ (P^) x 
with A ~ 7.5, but there is only a very weak correlation 
between V^/V to t and the packing probabilities. For ex- 
ample, the scatter in the data can vary by more than 
20 orders of magnitude! Thus, features of the basin ge- 
ometrical structure beyond the core region control the 




FIG. 7: (a) The weighted basin profile functions S^(Z) (for 
each of the Af s = 80 distinct MS packings for N = 6) sampled 
on hyper-spherical shells a distance I from MS packing n using 
method 2 with 6 = 1. (b) The scaled weighted basin profile 

function Sl(l) = [S£ (l)6r(k)e' l/e ] 1/{k ~ 1) plotted versus the 
scaled distance 1 = 1/0 for the same data in (a). The solid 
line has slope 1. 

MS packing probabilities for packings that are generated 
from dilute initial configurations. 

To begin to investigate the nature of the basin mor- 
phology beyond the core region, we characterize in detail 
the shapes of the weighted basin profile functions for each 
of the Af s MS packings for N = 6 in Fig. [7] (a). As found 
for the distribution of Voronoi volumes in dense granu- 
lar packings [l7l - tl9| . the form of S^(l) is described by a 
T-distribution 

where 9 = ((I 2 ) - (I) 2 ) /(I), k = (l)/0, and (I) = 
Jq 00 dllSn(l). The scaled weighted basin profile functions 

S r n (l) = [S'^(/)6T(fc)e- i / e ] 1 /( fe - 1 ) for all microstates col- 
lapse when plotted versus the scaled distance I — 1/6. 
The wider scatter at large I is caused by under sampling 
low probability configurations. We find similar quality 
for the collapse at larger N. 

We investigate the system size dependence of the av- 
erage weighted basin profile function 

Mb 

S r (l)=Y, P n S n( l ) ( 8 ) 
n=l 
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FIG. 8: The average weighted basin profile function S T (l) for 
several system sizes N = 2, 4, 6, 8, 10, 12, and 100 (from left 
to right) for damping parameter 6=1. The inset shows the 
parameters k (circles; left axis) and 1006 (squares; right axis) 
that describe fits of S T (l) to the F-distribution (Eq.O versus 
JV on a log-log scale. 



in Fig. |U over the range N = 2 to 100. S r (l) shifts to 
larger I with increasing N; the peak position k increases 
by a factor of 5 and scales roughly as \f~N over this range 
in N . The width 9 slightly narrows over the same range 
of N, scaling roughly as iV -1 / 4 . 

We also investigated the protocol dependence of the 
basin profile functions by varying the damping param- 
eter (b in Eq. [T]) used in the packing finder in method 
2. Lowering b decreases the rate at which energy is re- 
moved from the system and allows the system to explore 
larger regions of configuration space. In contrast, larger 
b increases the rate at which energy is removed from the 
system, and thus the initial configurations are typically 
closer to the final MS packings. In the top panel of Fig.Hl 
we plot the average weighted basin profile function ver- 
sus the damping parameter employed in method 2 over 
three orders of magnitude in b from 10 -2 to 10. We were 
able to saturate the b dependence of S r (l) for both large 
and small b, i.e. for b < 10~ 2 and b > 10, S T (l) is very 
weakly dependent on b. The two parameters k and 9 
that describe the shape of S r (l) exhibit two key features 
in the bottom panel of Fig. |9j 1) the peak of the dis- 
tribution (captured by k) and thus the lengthscales that 
determine the MS packing probabilities increase with de- 
creasing b and 2) the variance (relative to the average) 
depends weakly on b, but does possess a small peak near 
b = 10 _1 . We expect that the qualitative features of the b 
dependence will persist for larger system sizes. In future 
studies, we will predict the locations of the peaks in S T (l) 
at large and small b by calculating the distances between 
random points in configuration space and between ran- 
dom points and configurations related by particle label 
permutations, respectively. 




b 



FIG. 9: [top] The average weighted basin profile function 
S r (l) for N = 4 plotted over a wide range of the damping 
coefficients b employed in method 2. b ranges from 0.01 to 
10 from left to right, [bottom] The parameters k (circles; left 
axis) and 1006 (squares; right axis) that describe fits of S T (I) 
in (a) to the T-distribution (Eq. [7]) versus b. 



IV. CONCLUSIONS 

In this manuscript, we described and carried out a 
novel computational method for calculating the volume 
of the MS packing 'basins of attraction', which we de- 
fine as the collection of initial points in configuration 
space at zero packing fraction that map to a given MS 
packing by following a particular dynamics in the den- 
sity landscape. Note that our definition of the basin 
of attraction is protocol-dependent, and thus the basin 
volume will vary with the rate at which energy is dissi- 
pated, the compression rate, and other parameters. Us- 
ing dilute configurations as initial conditions and includ- 
ing variations in the basin volume with changes in the 
packing-generation protocol are crucial for understanding 
the protocol-dependent structural and mechanical prop- 
erties of granular media and other athermal particulate 
systems. 

Our computational studies of the basin volumes have 
uncovered three important results: 1) A small approxi- 
mately hyperspherical region of the basin of attraction 
with radius l c surrounds each MS packing, but the vol- 
ume of this region (relative to Vtot) is much smaller and 
only very weakly correlated with the MS packing prob- 
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FIG. 10: 113 snapshots of one of the N = 6 mechanically 
stable packings (shaded gray in Fig. 2} generated from inde- 
pendent random initial conditions. This packing contains one 
rattler particle (labeled 2) that can be positioned in the cav- 
ity on the left or right and at multiple positions in the right 
cavity. 

abilities in contrast to previous studies of jammed sys- 
tems 2) the probabilities of MS packings initialized 
with dilute configurations are instead controlled by fea- 
tures of the basins of attraction at lengthscales much 
further from the MS packing than the core region. In 
addition, the lengthscales that control the MS packing 
probabilities grow with increasing system size and de- 
creasing damping parameter b; and 3) The shape of the 
basin profile functions are well characterized by a T distri- 
bution, which suggests that we can construct a statistical 
mechanics- like theory to predict the shape of S r (l) and 
the MS packing probabilities. 

Acknowledgments This research was supported by the 
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Appendix A: Rattler Particles 

As shown in Figs. [3] and EJ MS packings contain rattler 
particles. Two of the Af s — 6 distinct MS packings for 
N = 4 and 24 of the Af s = 80 distinct MS packings for 
N = 6 contain rattler particles. For these small-TV sys- 
tems, the fraction of MS packings that contain rattlers 
is larger than the fraction of particles (roughly 5%-10%) 
that are rattlers in large MS packings [2fJ, but these re- 
sults suggest that the number of MS packings containing 
rattlers is extensive with Af s [21|. How do rattler parti- 
cles affect the calculation of the basins of attraction for 
MS packings? 

We find that the correspondence between the un- 
weighted fn(l) and weighted S„(l) basin profile functions 
breaks down for small / for MS packings that contain 
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FIG. 11: (color online) The basin profile function f^(l) on 
linear-log [top] and log-log [bottom] scales measured using 
methods 1 (circles) and 2 (solid lines) for the MS packing in 
Fig. [10] for 6 = 1. The vertical lines indicate l c . 



rattler particles. As shown in Fig. [TU1 for MS packings 
containing rattlers it is difficult to define uniquely the 
distance from the initial state to the final MS packing 
because the rattler particle can exist over a range of po- 
sitions for a given distinct MS packing. Further, the dif- 
ferent rattler locations may give widely varying contri- 
butions to the MS packing probability. 

In Fig. [TTJ we plot /£(£) calculated using methods 1 
(circles) and 2 (solid lines) for the MS packing depicted 
in FigfTTJl As described in Sec. lIIIl for method 1, we mea- 
sure the fraction of times the system returns to the initial 
MS packing after a perturbation of size I, which is largely 
unaffected by the presence of rattlers. For method 2, we 
measure the normalized distribution of distances between 
the initial configurations and the final MS packings. For 
large I, method 2 is also largely unaffected by the pres- 
ence of rattlers. However when the initial configuration 
and final MS packing are close together (i.e. small I), the 
fact that the rattler is not always in the same position in 
the final MS packing leads to a significant error in mea- 
suring I and hence /^(Oi as shown in Fig [IT] For our 
measurements of basin volumes at small I, such as V£ in 
Fig. [51 we show results using method 1. Our main results 
are insensitive to the presence of rattler particles because 
MS packing probabilities (generated from initially dilute 
configurations) are determined by features of f„(l) at 
large I. 
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